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Abstract 

A Lattice Boltzmann model is considered in which the speed of sound 
can be varied independently of the other parameters. The range over which 
the speed of sound can be varied is investigated and good agreement is 
found between simulations and theory. The onset of non-linear effects due 
to variations in the speed of sound is also investigated and good agreement 
is again found with theory. It is also shown that the fluid viscosity is not 
altered by changing the speed of sound. 

1 Introduction 

The lattice Boltzmann model (LBM) is a numerical technique for fluid simula- 
tion which has become increasingly popular in recent years [H [21 El S] ■ The LBM 
originates from the lattice gas model (LGM) O El [71 [8] were fluid particles are 
constrained to move on a regular lattice such that their collisions conserve mass 
and momentum. The particles are further constrained to move with unit veloc- 
ity and with a maximum occupancy of one particle per grid direction per grid 
site. The evolution of the LBM from the lattice gas model involved a number of 
key developments. The Boolean particle number (1 if a particle is present and 
otherwise) was replaced by a real number, later recognized as the distribution 
function, representing an ensemble average of the particle occupation [9J; this 
removed the noise associated with the relatively small number of fluid particles 
represented in the LGM. Obtaining an ensemble average of the particle collisions 
is not straightforward but the process was simplified to depend only on the link 
direction [10] and the isotropy of the model [IT]. Finally the particle-particle 
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collisions were replaced by considering the distribution functions relaxing toward 
a Maxwell-Boltzmann equilibrium distribution |12j . 



The LBM is referred to as an incompressible technique because the LBM 
scheme can be shown to satisfy the incompressible Navier-Stokes equation in the 
limit that the density, p, does not vary in space or time. However, when apply- 
ing the LBM, there is no restriction on the density to remain constant and in 
many applications the density will vary in space and/or time. In practice the 
'incompressible' nature of the LBM is interpreted as requiring that any density 
variations are small or that the density varies only slowly in space and time. In 
this limit it is possible to simulate phenomena such as acoustic waves where a 
density variation is required, provided the variation is small [131 [HI fT5] . 

In the LBM, pressure is defined through the equation of state, see for example 
[2], p = c^p, where c is the speed of sound which takes a fixed value in any sim- 
ulation, dependent only on external factors such as the shape of the simulation 
grid. Thus we see that the speed of sound controls the relationship between the 
density and the pressure and therefore the compressibility of the fluid. To model 
fluids with different compressibility we require to be able to change the speed of 
sound in the simulation. Here we consider a model with a variable speed of sound 
and investigate the non-linear aspects of sound wave propagation. 



2 The lattice Boltzmann model 

In this section we briefly consider the standard LBM on a square grid in two 
dimensions. In the LBM the distribution functions, fi{x,t), for i= 0-8 evolve on 
a regular grid according to the Boltzmann equation [16], as 



where = (cos[7r(i — l)/2], sin[7r(i — l)/2]) for i = 1-4, and = -\/2(cos[7r(2 — 
9/2)/2], sin[7r(i — 9/2)/2]) for i = 5-8 are link vectors on the grid and Bq = (0,0). 
The left hand side of equation ([1]) represents streaming of the distribution function 
from one site to a neighbouring site. The right hand side is the Bhatnagar, Gross 
and Krook (BGK) collision operator [171 [121 [iH]- The equilibrium distribution 
function is given by |16] 



fi{x + ei,t + l) - fi{x,t) 
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fi = Wip l + 3ei-u + -{ci ■ uf 
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where the fluid density, p, and velocity, u are found from the distribution function 
as 

8 8 

P = ^fi and pu = J2 ft^i (3) 

i=Q 1=0 

and where Wq = 4/9, Wi = 1/9 for i=l-4, and Wi = 1/36 for i = 5-8. The 
relaxation time, r, determines the rate at which the distribution functions relax 
to their equilibrium values; it determines the fluid viscosity, z/, as 

Equation ([T]) can be modified by the addition of an extra source term. This 
approach has been used to implement a body force [19] and also to simulate 
non-ideal equations of state such as phase separation [20] . 



3 A lattice Boltzmann Model with a Variable 
Speed of Sound 

One method for achieving a variable speed of sound was proposed by Alexander et 
al. [5T] who considered a model on a hexagonal grid with an altered equilibrium 
distribution function in which the ratio of 'rest particles' (/o) to 'moving particles' 
(/j, ^ 7^ 0) can be altered. Simulations performed using this model [21j showed 
that the speed of sound can be varied between and approximately 0.65 with 
good agreement between theory and simulation for speed less than approximately 
0.4. Above 0.4 there is some deviation between theory and simulation. In this 
model the viscosity was also found to be a function of the variable speed of sound. 

An alternative approach was proposed by Yu and Zhao [22j. They intro- 
duced an attractive force which produced a variable speed of sound which was 
dependent on the amplitude of the introduced force. The model was verified by 
measuring the Doppler shift and the Mach cone for Mach numbers less than and 
greater then unity respectively. 



Here we consider the following LBM: 

f,{x + ei, t + 1) - f,{x, t) = -!(/, - 7J + 3wiaVp ■ e^, (5) 

r 

where f ^ and Wi are defined as previously and the additional term on the right- 
hand side of equation ([5]) represents a body force aVp, see for example [19]. 
This is effectively the same as the model of Yu and Zhao [22], except that the 
forcing term is expressed explicitly as proportional to the density gradient and 
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the amplitude term, a is not restricted to be positive. Performing a Chapman- 
Enskog expansion [8] we can write the Navier-Stokes equation in the form 



a 



da [cl -a) p + ud(}{daUp + dfsUa), 
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where 
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as before, and we have assumed that the derivatives of p can be neglected except 
where they appear in the pressure term. This is the Navier-Stokes equation for 
a fluid which has pressure p = c^p, where Ce is the effective speed of sound which 
is given by 



Thus by introducing an additional term to the Boltzmann equation which acts as 
a body force proportional to the density gradient, we have included a force which 
increases or decreases (depending on the sign of a) the pressure forcing term in 
the Navier-Stokes equation, and hence the speed of sound in the model. 

A number of models, based on the cellular automata (CA) approach of the 
LGM, have also been applied to simulate acoustic waves as well as Burgers' 
equation. Mora [23] introduced the phononic lattice solid model which obeyed a 
Boltzmann equation and satisfied the acoustic wave equation in the macroscopic 
limit. The Boltzmann equation was solved using a finite-difference scheme. This 
model was further developed [21] in the phononic lattice solid by interpolation 
model in which the particle number densities move along the lattice links in the 
same manner as the distribution functions in the LGM and LBM. Unlike the 
lattice gas particles they travel at different speeds and so, in general, do not 
arrive at a grid site after each time step. The number densities are therefore 
interpolated to find the value at each site. 

CA models for Burgers' equation have also been considered. The model of 
Boghosian and Levermore [25] was based on the LGM approach and was shown 
to follow the solution of Burgers' equation. The convergence of the model was 
also established [26]. Following the development of the LBM from the LGM 
[Ol [10, [H] , Elton [27j considered a CA model for Burgers' equation based on the 
mean occupation number rather than discrete particles. The model compared 
favourably with a finite difference solution of Burgers' equation. A quantum lat- 
tice gas model [28] and an intrinsically stable entropic model [29] have also been 
proposed recently. Each of these models provides a simulation method which 
satisfies Burgers' equation. This is different to the model investigated here which 
satisfies the Navier-Stokes equation and is also capable of simulating non-linear 
acoustic waves whose behaviour, for the case a = 0, has been shown to be in 
good agreement with Burgers' equation [T3]. This means that it is capable of 
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simulating the interaction between acoustic waves and fluid flow, see for example 
[I5l[30l[3ll[32]. 



4 Numerical Simulations 

The variation in the speed of sound was investigated in the following simulations 
of plane acoustic wave propagating in an unbound media. This if effectively a 
one dimensional problem so a gird was used with only 4 points in the ^/-direction. 
The grid had L points in the x-direction and was initialised with zero velocity 
and unit density. Periodic boundary conditions were applied at y = and y = 
4. At a; = a boundary condition was applied in which the density was varied 
in a sinusoidal manner with a period of 500 time-steps, the velocity at this grid 
boundary was maintained at zero. At y = L a boundary condition of p = 1 
and u = was applied. L was selected to be large enough so that the density 
disturbance does not reach x = L during the measurement window. From these 
simulations the speed of sound was found in two ways: 1) from the time taken for 
the wave to pass between two points at a known separation; and 2) since we are 
dealing with a low amplitude plane harmonic wave, Ce can be derived from [33] 
u{x*,t*)/ce = [p{x*,t*) - po]/po, where u{x*,t*) and p{x*,t*) are the fluid ve- 
locity and density respectively, at some position x* and time t*; here we selected 
X* and t* such that they correspond to a pressure and velocity maximum. The 
results are shown in figure [T] where there is excellent agreement between the mea- 
sured values of Ce and the theoretical values given by equation over a range of 
a. At « = 1/3 we have Ce = giving an upper limit for a. Now a = —2/3 gives 
Ce = 1 which would appear to be an upper limit since this is the speed with which 
the distribution functions propagate and so is the maximum speed information 
about the density, pressure and velocity can be transmitted through the fluid. 
Initial simulations were performed using a simple forward difference scheme to 
calculate the pressure gradient. In this case the model became unstable as Ce 
approached unit, as would be expected. When a central difference technique was 
applied to calculate the density gradient the model was found to remain stable 
as Ce approaches unity. Indeed it was found to remain stable for higher values 
of Ce (a < -2/3) and the simulated density waves were observed to propagate 
with Ce > 1. The transfer of information at a speed greater than unity can be 
understood with the introduction of a central difference scheme to calculate the 
'pressure' forcing term. Consider a density distribution where p{x) = po every- 
where at t=0, except for p{xo) = po + 6p. Then at x = xo + 1 the density gradient 
(calculated using a central difference scheme) will be non-zero. In the case of a < 
this will give a forcing term that will cause the value of the distribution function 
travelling in the positive direction at Xq + 1 to increase. This information (about 
a density increase at xq) will then be passed on to xo-l-2 at t=l by the streaming 
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action. As seen in figure [T] this is a limited effect; however, it suggests that LB 
schemes in which the speed of sound is significantly greater than unit may be 
possible if a longer range action is applied. The largest value of Ce which was 
observed to remain stable was Ce = 1.125. The validity of simulations with Cg > 
1 has not been established and requires further investigation. The remainder of 
the simulations presented here have have Ce < 1. 

To further investigate the effectiveness of the model two progressive sound 
waves with A = 4000 were initialised with the same velocity amplitude: u{x,t = 
0) = Mosin(fcx), v{x,t = 0) = and density profile p*^*^(x,t = 0) = po\l + 
u^/c^e^) sm{kx), for z = 1, 2; where c^*^ is found from a'-*-' with a^^^ = 0.2933 and 
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-0.47667. This gives = 0.2 and 4^) = 0.9 such that c^^)/^^) = 4.5. 
The ambient density pg ^ is also free to be varied, however, in the simulations we 
used p'q^ = p'q^ = 1 and uq = 0.002. For these two waves we can calculate the 
shock development distance given by 

where M'^*-' = -u/c^*^ is the Mach number, k = 2n/\is the wavenumber and e = 1. 
This is the distance over which an initially sinusoidal wave in an inviscid fluid 
will develop into a discontinuous wave with a sawtooth shape. Here we have 
^(1) ~ iQ\ Q^Yid ^(2) so we expect wave (1) (in the fluid with the lower 

speed of sound) to exhibit stronger non linear effects than wave (2). 



Following [T^ periodic boundary conditions were applied on all grid bound- 
aries and the waves were allowed to propagate. Figure [2] shows the normalised 
velocities u' = u/uq; plotted as a function of = uit — x/ce) after the waves 
has propagated distances of 5A, lOA and 15A. Here u is the angular frequency 
of the acoustic wave. It can clearly be seen that wave (1) exhibits strong non- 
linear behaviour over the first twenty wavelengths. Also shown in figure [2] is the 
numerical solution of Burgers' equation [14] which describes the wave form in a 
viscous fluid in terms of a truncated sum of N harmonics. Good agreement is 
observed for x = 5, 10 and 15 A. For x = 20A a suitable value of N could not be 
found to give a smooth profile. Even at large N oscillations were observed on the 
wave form around = and = ±1. For this reason the solution of Burgers' 
equation is only shown away from these points for a; = 20 where there is again 
good agreement with the Boltzmann simulations. The LBM simulation of wave 
(2) at X = 45A(= 10(4^-*/^^^^) A) is also shown in figure [2] for comparison. The 
wave form is similar to that for wave (1) at x = lOA as would be expected from 
equation (Q since both waves were run with a relatively low viscosity (z/ = 0.01). 
A comparison with wave (1) is given in figure [3] which shows the amplitudes, a„, 
of the fundamental {n = 1) and higher (n > 1) harmonics. The results in figure 
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[3] were obtained by performing a fast Fourier transform on the velocity between 
X = 14A and x > 20A where the smallest number of point which were an exact 
multiple of 2 were used. 

Figure H] shows the long term behaviour of wave (2) for three different fluid 
viscosities given by z/ = 1/30, 1/2 and 7/6 , the value of A was 400. Relatively 
large values of r were selected to ensure that the wave was damped rapidly. The 
results are in good agreement with the exponential decay rate predicted by linear 
theory, see [Ml ESj [13] , which is also shown. 

The long term behaviour of wave (2) was also considered in fluids with a lower 
viscosity. The decay rate per unit time is independent of the speed of sound and 
depends only on the wavelength of the wave and the viscosity of the fluid. This 
is shown in figure O which depicts wave (2) for a range of a values after around 
100 periods for the wave with a = —0.65. It is clear from figure [S] that the 
damping rate is the same for each value of a; including a = which corresponds 
to the standard LBM with no additional body force. The viscosity of the fluid 
was obtained from the simulations as 



and is shown in figure [6]for a = —0.6 for values of r between 0.501 and 0.9. Fig- 
ure [6] shows good agreement between the measure viscosity and the theoretical 
value at large viscosities. At lower viscosities their is some deviation. For a = 
-0.6 the largest difference observed between the simulated wave amplitude after 
100 periods and the theoretical value was less then 4 %. It is well known that 
the LBM becomes unstable as r ^ 0.5 with noise being introduced into the sim- 
ulation which eventually becomes unstable. This has been observed elsewhere, 
for example p3] where a filter was introduced to reduce the noise. Here a filter 
was not required and the introduction of the body force term was not observed 
to significant alter the onset of instabilities. The good agreement between the 
simulations and the theory indicates that the introduction of a variable speed of 
sound using the method proposed here does not effect the viscosity of the fluid, 
as predicted by equations and ([7]). 



5 Conclusions 

A new lattice Boltzmann model has been devised in which the speed of sound can 
be varied. This has been demonstrated in a number of simulations in which the 
speed at which a disturbance propagates in a fluid was found to agree with the 
theoretical speed of sound. Further, it has also been shown that the non-linearity 
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of a pressure wave depends not only on its amplitude but also on the speed of 
sound in the medium through which it is propagating, as would be expected, 
and good agreement was found with the numerical solution of Burgers' equation. 
Compared to a previous model [21], the present model has a larger range over 
which the speed of sound can be varied and shows better agreement with theory 
over the full range. Additionally the method by which the variable speed of sound 
is introduced does not change the viscosity of the fluid. The speed of sound is 
varied in the new model by a free parameter, a, which was kept constant in 
space and time during the simulations presented here. This is not, however, a 
fundamental requirement and it is envisaged that this model could be applied 
such that the speed of sound varies within the simulation. For example, in a 
simulation of a liquid-gas or a binary fluid mixture, the speed of sound could be 
varied as a function of the fluid species. This could be done in much the same was 
as a species dependent relaxation time has been used to simulate binary fluids 
with different viscosities, see for example [SB]. Alternatively, in a single phase 
fluid simulation, the speed of sound could be varied in a pre-determined manner 
to account for an external influence; for example a temperature gradient or an 
acoustical lens. The model can also be applied to simulate fluids with different 
compressabilities which extends the range of application of the LBM. 
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Figure 1: The effective speed of sound, Ce as a function of the forcing magnitude 
a. 
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Figure 2: The nonlinear development of the sound wave for Ce = 0.2 over the 
first twenty wavelengths. Also shown for comparison is the numerical solution of 
Burgers' equation (+) and the LBM simulation of a sound wave with Cg = 0.9 at 
X = 45A (diamonds). 
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Figure 3: Amplitude of the spectral components for wave (1) and wave (2) cal- 
culated over just over six wavelengths starting at x — 14A. 
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Figure 4: The propagation of wave (2) for three selected values of r. Also shown 
is the e^ 




Figure 5: Propagation of wavc(2) for r = 0.52 and selected values of a for 
97y(-o.65) < ^ < 100T(-°-65), where T^-^-^^) is the period of the wave when 
a = -0.65. 
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Figure 6: The fluid viscosity calculated from the decay of a wave after 100 periods 
with a = -0.6, as a function of the relaxation parameter. The solid line represents 
the theoretical expression. 
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